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METHOD FOR PACKING AND FILTERING GEOPHYSICAL EVENTS TO 
RETRIEVE THEREFROM DATA ON THE TYPE OF SUBSOIL 



The present invention relates to a method for packing and filtering geophysical data 
in order to retrieve therefi'om information on the nature of the subsoil. 

This methodology is applied to packing and filtering sets of multi-donGuun 
geophysical data, i.e. either pre-stack seismic data fi'om a 3D seismic survey grouped 
into a certain number of data cubes corresponding each, for example, to the same offset 
(seismic trace cubes referred to as iso-ofifset cubes) or to the same angle of incidence 
(referred to as iso-incidence angle cubes), or seismic attribute cubes obtained by seismic 
amplitude processing such as, for example, elastic properties cubes firom a joint 
stratigraphic inversion of pre-stack seismic data, these cubes being obtained by 
implementing the method described for example in patent FR-2,800,473 (EP- 1,096,270) 
filed by the applicant. It then consists in shaping the pre-stack seismic information 
contained in seismic trace portions formed by exploration of an underground zone, by 
means of new representative attributes which are thereafter used for geologic 
interpretation or interpretation of an underground reservoir. It can also apply to packing 
and filtering of pre-stack or post-stack 3D seismic data fi-om repeated surveys (4D 
seismic method) in order to interpret the fluid content variations of a petroleum 
reservoir under production. This method is also applicable to packing of log data of 
several physical natures by forming, fi-om initial measurements, new log attributes in 
order to improve interpretation of the log data in terms of lithologic and petrophysical 
characteristics of the reservoir. 



FIELD OF THE INVENTION 
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BACKGROUND OF THE INVENTION 

Geophysical measurements are conventionally used to provide complementary 
information, in relation to direct observations in wells, on the geologic properties 
variations of underground formations. Seismic interpretation is for example based on 
5 the analysis of attributes in order to retrieve these geologic properties from the seismic 
information, in a zone of interest of an underground formation (reservoir zone for 
example). Another example is log interpretation where the various measurements are 
simultaneously analysed to deduce therefrom the petrophysical characteristics of the 
underground formation studied. 

10 A standard approach for seismic or log attributes analysis involves the use of 

statistical pattern recognition and estimation methods in order to calibrate the pertinent 
attributes with the available information obtained in various wellbores. 

For example, patent FR-A-2,768,818 (US-6,05 1,651) ffled by the applicant 
describes a method allowing to detect the potential classes in a population of events 
1 S related to the physical properties of a complex medium such as the subsoil, located from 
data obtained by subsoil exploration or by in-situ measurements (events of geologic, 
geophysical nature, events linked with fluid production data, etc.), these events being 
located by points in a multi-variable space and defined by an attributes vector. 

Another example is given by patent FR-01/05,675 filed by the applicant, which 
20 describes a method for facilitating monitoring with time of the evolution of an 
undergroimd zone by compared analysis of a certmn number n of seismic records 
successively obtained after n successive 3D seismic stirveys (4D seismic method), 
wherein a pattern recognition technique applied to the whole of the seismic events of 
several surveys, considered and analysed simultaneously, is used. 
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Patent EP-A-671,017 (US-5,63 8,269) filed by the applicant can also be mentioned, 
which describes a method allowing to establish a connection between geologic data 
obtained by coring or logging in wells and seismic data obtained by means of a seismic 
exploration survey, which is based on a statistical calibration technique with statistical 
5 calibration obtained by bringing together the local geologic data (measured in wells) 
and seismic attributes read on seismic traces obtained in the immiediate vicinity of each 
well. 

The seismic attributes that are subjected to such an intapretation are calculated in 
the post-stack or in the pre-stack domain. 

10 In the post-stack domain, the conventional attributes are calculated fi*om the 

amplitudes at the reservoir level, or firom the impedance P estimated by means of a 
stratigraphic inversion. 

In the pre-stack domain, the number of pertinent attributes can increase 
considerably. In this case, the attributes are either the amplitudes at different ofi&ets or 
15 incidence angles, or parameters obtained from a pre-stack joint stratigraphic inversion 
(impedances P and S, density, product of the Lame parameters with the density, etc.). 

Concerning the interpretation of logs, the situation is similar, a very large number of 
measurements being available, hence a large number of attributes allowing to 
characterize the logs, these attributes being often organized in families of equal physical 
20 significance (for example fiunily of the resistivity attributes, family of the attributes 
associated with the natural or induced radioactivity of the formations, etc.). 

A method referred to as principal components analysis (PC A) method, well known 
to the man skilled in the art, can be applied to these attributes in order to analyse the 
relations existing between them and to reduce the number of significant attributes. The 
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principal components extracted from the PCA are the new attributes : they define an 
orthogonal or non-orthogonal (rotation) vector base and are linear combinations of 
initial variables. They can be used in statistical pattan recognition algorithms and 
correspond to a multivariate filtering of the initial seismic or log information. Various 
5 appUcations of the PCA method are for example described in : 

- Dumay, J., Foumier, F., 1988, 'Multivariate statistical analyses applied to seismic 
facies recognition". Geophysics, 53, 1 151-1159, 

- Hagen, D., C, 1982, "The application of principal components analysis to seismic data 
sets", GeoexpL, 20, 93-111, or 

10 - patent FR-A-2,772,138 (US-6,345,108) filed by the applicant. 

However, PCA does not take account of the notion of attributes group and, 
consequently, it is ofi:en difiQcult to ^ve a physical interpretation of the principal 
components, or to relate them clearly to the initial attributes, especially if they are 
organized in groups of equal physical significance. 

15 The analysis referred to as canonical analysis is also a statistical method well 

known to the man skilled in the art, which allows to establish the relations that may 
exist between two sets of variables in order to know if these two sets describe the same 
properties. This method is for example described by : 

- Hotelling, H., 1936, "Relations between two sets of variables", Biometrika, 28, 321- 
20 377. 

An example of application of canonical analysis is for example described by : 

- Foumier, F., and Derain, J.F., 1995, "A statistical methodology for deriving reservoir 
properties from seismic data". Geophysics, 60, 1437-1450. 
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This method is limited to the study of two sets of variables and does therefore 
generally not apply to the multicube seismic information or to the miiltidomain log 
information. Furthermore, it defines synthetic variables in each subspace associated 
with the two sets, and not a single vector base allowing to pack the whole of the initial 
5 variables and to describe the relation between the two sets. 

There are several methods allowing to generalize the analysis referred to as 
canonical by extending the analysis to more than two sets of variables. Different aspects 
of the prior art in the field considered are described, for example, in the following 
publications : 

10 - Horst, P., 1961, '"Relations among m sets of measures", Psychometrika, 26, No.2, 129- 
149, 

- Caroll, J.D., 1968, "A generalisation of canonical correlation analysis to three or more 
sets of variables". Proa 76**^ Conv. Amer. Psych. Ass., 

- Kettering, J.R., 1971, "Canonical analysis of several sets of variables", 58, 3, 433- 
15 450, or 

- Saporta, G., 1990, 'Trobabilites, analyse des doimees et statistiques", Technip, Paris. 

These methods however have certain limits. In particular, they do not allow to 
describe each set separately, the synthetic variables cannot be readily related to the 
various sets and, consequently, be physically interpreted. Furthermore, only the global 
20 inertia part they represent can be known, and not the inertia part they represent for each 
set. It is therefore not possible to filter each set independently. 
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Another known analysis technique referred to as generalized principal components 
analysis (GPCA) allows to conopare diflferent sets of variables while describing each 
one. It is implemented on data of economic nature for example in : 

- Casin, Ph., 2001, "A generalization of principal component analysis to K sets of 
5 variables". Computational Statistics & Data Analysis, 35, 417-428. 

SUMMARY OF THE INVENTION 

The method according to the invention provides an approach wherein, by 
implementing a multivariable statistical technique applied to several sets of seismic 
records or logs available, obtained by e?q)loration of an underground zone, the 
10 geophysical events readable on these sets of measurements are packed and filtered. The 
method comprises : 

- distributing these geophysical events into a certain number p of families of 
geophysical attributes having each a particular physical significance, 

- forming, by combination of said geophysical attributes, synthetic variables in much 
IS fewer number than the number of said attributes, which are obtained by constmction of 

an orthogonal vector base in each one of the p analysis sets consisting of the data of 
each one of said &milies, wherefi-om is derived the formation of an orthonormal vector 
base describing the p analysis sets, 

- using this orthonormal vector base for filtering and describing said geophysical events. 

20 Formation of the synthetic variables is obtained for example by combining a 

principal components analysis technique or its d^ved techniques based on 
diagonalization of a norm (seeking eigenvectors ) and a simple regression technique. 
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According to an implementation mode, the synthetic variables can be subjected to a 
rotation so as to improve the correlations between initial variables, i.e. the geophysical 
attributes, and synthetic variables. 

The synthetic variables can be formed from families comprising for example one or 
5 more seismic trace cubes obtained after a 3D seismic survey, each one corresponding to 
the same offset. 

Synthetic variables can also be formed from families consisting of seismic trace 
cubes obtained after a 3D seismic survey, each one corresponding to the same angle of 
incidence. 

10 Synthetic variables can also be formed from families consisting of seismic attribute 

cubes formed by processing the seismic trace amplitudes obtained at the end of a 3D 
seismic survey. 

Synthetic variables can also be formed from families consisting of seismic trace 
cubes obtained after successive seismic surveys in the zone. 

15 It is also possible to form synthetic variables from families consisting of seismic 

attribute cubes formed by processing the seismic trace amplitudes obtained after 
successive seismic surveys in the zone. 

Synthetic variables can also be formed from families consisting of measurements of 
log attributes of different physical nature, such as the family of resistivity 
20 measurements, the family of measurements related to the formation radioactivity, the 
family of measurements related to the formation porosity, etc. 

The method is particularly advantageous in that it allows, when applied to seismic 
or log data : 



- to define new geophysical attributes, linear combinations of the initial geophysical 
attributes, summarizing the initial seismic or log data and allowing to filter this data, 
while respecting the multidomain structure of the data, 

r to show if all the families of initial geophysical attributes, or only some of them, 
5 provide redundant information, and to measure the part of this redundant information, 

- to give a physical interpretation to the new geophysical attributes. 

The new geophysical attributes thus formed can then be analysed in order to 
retrieve the geologic properties of the seismic or log data. Two classes of methods can 
be conventionally used therefore : a first one using statistical estimation techniques and 
10 providing quantitative interpretation in terms of geologic properties, and a second one 
using pattern recognition techniques to analyse the seismic facies or the log facies 
(electrofacies) and provide qualitative information. 

Examples of these two types of methods are respectively implemented in the 
aforementioned patents EP-A-671,017 (US-5,63 8,269) and FR-A-2,768,818 (US- 
15 6,051,651) filed by the applicant. 

BRIEF DESCRIPTION OF THE FIGURES 

- Figures lA to IC show the organization of the seismic records in p sets : each set 
groups in a cube all of the trace portions corresponding to ah identical physical 
criterion, 

20 - Figures 2 A to 2C respectively show the maps of the first three new seismic attributes 
Z(^>,Z<^>andZ<^>, 
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- Figures 3 A to 3D respectively show the first true temporal plane (3 A) of the analysis 
window and the same plane reconstituted/filtered with 1 (3B), 3 (3C) or 8 (3D) of the 
new seismic attributes for the cube of impedance P, 

- Figures 4A to 4D respectively show the first trae temporal plane (4A) of the analysis 
5 window and the same plane reconstituted/filtered with 1 (4B), 3 (4C) or 8 (4D) of the 

new seismic attributes for the cube of impedance S, 

- Figures 5 A to 5D respectively show the first tme temporal plane (SA) and the same 
plane reconstituted/filtered with 1 (5B), 3 (5C) or 8 (5D) of the new seismic attributes 
for the cube of density p, 

10 - Figures 6 respectively show the vertical sections, passing through well W2, tme (6A), 
reconstituted/filtered with the first 8 new attributes (6B), and the difiTerences (residues) 
(6C) between these two sections for the cube of impedance P, 

- Figures 7 respectively show the vertical sections, passing through well W2, trae (7A), 
reconstituted/filtered with the first 8 new attributes (7B), and the differences (residues) 

1 5 (7C) between these two sections for the cube of impedance S, 

- Figure 8 shows the location of the learning traces used to calibrate the classification 
function within the context of the supervised pattern recognition of the firat 8 new 
seismic attributes, 

" Figures 9A, 9B respectively show the assignments of the seismic events and the most 
20 reliable assignments of these events within the context of the supervised pattern 
recognition of the first 8 new seismic attributes, and 



10 



- Figures lOA, lOB respectively show the assignments of the seismic events and the 
assignments of these seismic events within the context of the non-supervised pattern 
recognition of the first 8 new seismic attributes. 

DETAILED DESCRIPTION 

Several sets of geophysical records obtained by seismic e?q)loration or logging in an 
underground zone studied are available. The geophysical events to be packed/filtered 
are characterized by geophysical attributes often organized in sets or families of 
particular physical significance (Figures lA to IC). ITiese attributes can be the 
succession of the amplitude values of seismic trace portions or of elastic parameters 
fi-om a pre-stack joint elastic inversion, obtained for example by means of the method 
described in the aforementioned patent FR-2,800,473 (EP- 1,096,270) along the time 
window studied. These attributes can be, for a given depth, the measurements of a 
single set of logs obtained along a well, such as for example the values of the different 
resistivity tools, the values of the different tools measuring the porosity, etc. 

1 - Notations 

Let X be a matrix consisting of p juxtaposed tables Xi, . . . . . . : 

X= {Xi ... ... Xp} , Xj = pCj,i ... Xj^jg] 

where is a table consisting of mj variables Xjx of n individuals (in the statistical sense 
of the term). Matrix X groups together all the geophysical events available. These 
events are first organized in p sets or families Xj (Figure lA). Each family Xj 
corresponds to a particular physical property. For analysis of the 3D seismic survey, this 
family has the form of a cube of values corresponding to this property, located by their 
spatial position (Ime, cdp, time/depth) (Figures IB and IC). These cubes can for 



example consist of a set of amplitude values of n traces or trace portions associated with 
a common offset or a common angle or the spatial distribution of the values of a 
common attribute derived from these amplitudes such as, for example, the acoustic 
impedance in P or S waves, the density or Lame parameters, in a selected time window 
of size mj. In this example, each variable Xj^ groups together the values of the seismic 
attribute considered (seismic amplitude or elastic property) for a temporal sample within 
the analysis window. In this example, a cube comprises as many individuals as traces 
available. For analysis of the logs, this family has the form of a sequence of mj series of 
values, located by their position along the well (time or depth). This sequence can 
consist, for example, of all the logs of the tools of same physical nature, such as the 
various tools measuring the formation porosity, or the tools measuring various 
resistivity types, on a set of n time or depth points along one or more wells at the level 
of the undergroimd formation studied. In this example, each variable groups 
together the values of the log attribute considered (resistivity, porosity, radioactivity,. . .) 
for a time or depth point within the analysis window. In this example, a sequence 
comprises as many individuals as time or depth points available along the well(s) 
analysed. 

The mj variables of a table define a space denoted by Sj and are assumed to 
be centred. Let Dn be the diagonal matrix of the wdghts of the individuals. 

Let rij be an orthogonal projection operator on space Sj defined by : 
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2 - Calculation of the new attribute 

The purpose of Generalized Principal Components Analysis (GPCA) is both to 
describe the proximities between the p families of variables and to describe separately 
each set by means of synthetic variables that are linear combinations of the initial 
variables. Then, what is common to all the sets or families, or to some of them only, and 
consequently what is different can be shown. GPCA allows to reduce the number of 
variables used to describe a data set and it allows multivariate filtering while keeping 
the multidomain structure of the initial data. To reach this objective, the first synthetic 
variable calculated, called Z^^\ has to be such that : 

1. The p orthogonal projections Zj^*^ of iP^ on spaces Sj have to be as close as 
possible to one another; 

2. In each group j, the dispersion of the n individuals projected on Zj^^^ is maximum. 
This is mathematically translated by the following conditions : 

2.iMxgCov^(Zf,X,,,) 

where is the muhiple correlation coeflBcient and Zj^^^ is the projection of Z^^^ on 
space Sj normalized at 1. 

A new synthetic variable ZP^ is then sought according to the same criteria. An 
orthogonality condition is added between the projections of variables Z^^ of the same 
set. The procedure stops when an orthogonal base is found in each space Sj. The number 
of variables Z^^ does not exceed the largest dimension of spaces Sj. 
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Rotation methods well known to the man skilled in the art can be implemented in 
each space thus found in order to optimize the correlation coefficients between the new 
variables and the initial variables. 

The first variable tP^ is the normalized eigenvector of covariance matrix XJ^Dn 
corresponding to the greatest eigenvalue. Z^^^ is then projected on the p spaces Sj and the 
regression of the mj variables by Zf^^ is calculated in each space. The matrix of the p 
juxtaposed residues Xj^^^ defines the new space studied. This new space corresponds to 
the initial space from which the influence of Zf^^ has been removed in each set Sj. The 
mj variables of table Xj^^^ define a subspace Sj^^^ of Sj. The second component ZP is the 
normalized eigenvector of covariance matrix X^^^X^^^Dn of the residues corresponding 
to the greatest eigenvalue. ZP^ is projected on the p spaces Sj^^l The regression of the mj 
variables by Zj^^^ is calculated in each space Sj^^\ The regression carried out in each 
space Sj ensures that the projection of ZP^ on space Sj^^^ is orthogonal to the projection 
ofZ^^^ on space Sj. 

Calculation of the eigenvectors is carried out by means of any method known to the 
man skilled in the art. 

In practice, matrix XXton being of order n, where n can be very great, the 
normalized eigenvector V*^ of matrix X*^DnX^^ associated with the greatest eigenvalue 
is calculated in a stage k. Vector Z^^ is then simply calculated by linear transformation 
Z*> - x^H^) and normalized. 

Each set can be reconstituted with a small number of attributes. The initial variables 
of set j can be expressed as the sum of a linear combination of variables Zj^^ and of a 
residue. If all the variables Zj^^ of a set are used, the residue is zero. The coefficients of 
the linear combination are the coefficients calculated after each regression in spaces Sj. 
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One of the properties of the method according to the invention is that variables Z*^ 
are orthogonal two by two. 

For interpretation of the synthetic variables calculated, we then calculate for each 
set j: 

1. The square of the correlation between variable Z*^ and variable Zj^^ of set j, 
indicating the degree of proximity between these two variables; 

2. the inertia percentage of set j explained by variable Zj^^ . 

The correlations between variables ZP^^ and initial variables Xjjc are also calculated 
to show possible relations between the Xj^c. The correlations between variables Zj^^ and 
JQjc of the same set j, or the correlations between variables Zj^^ for the same k can also 
be calculated. 

The new synthetic variables thus created constitute the new geophysical attributes 
sought. 

3 - Application example 

The method according to the invention was applied to three cubes of elastic 
properties (impedance P denoted by Ip, impedance S denoted by Is and density denoted 
by p) obtained after pre-stack joint elastic inversion, by means of the method described 
in the aforementioned patent ER-2,800,473, of seismic records grouped in cubes of 
common incidence angles. Each elastic property generates a set of seismic attributes or 
variables. The analysis time window consists of 19 signal samples taken from seismic 
records. We consider the values of the elastic property for each temporal sample within 
the analysis window to be variables. Each set thus consists of 19 variables made up of 
as many individuals as there are traces available. Since there are three variable sets, we 
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have in total 57 variables. The three parameter cubes have been flattened in relation to a 
reference surface. 

The synthetic variables Z^^ obtained summarize the respective information of the 
three elastic properties, while taking into account the relations between them. 

5. 7 - Analysis of the relations between impedances P, S and density p 

Table 1 hereunder gives the variance proportion in percentage explained by the 
projections of the first 8 synthetic variables Z^^^ to Z^^^ in the three sets. It shows that 8 
variables are sufficient to summarize 87.4 % of the variance of impedance P, 94.4 % of 
the variance of impedance S and 88.4 % of the variance of the density. 
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Table 2 hereafter gives the square of the correlation between synthetic variables tP^ 
and their projections on the three sets. This value indicates the degree of proximity 
between the variable and the set. The closer the value is to 1, the closer the variable is to 
its projection and therefore to the set. 
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For example, 2^*^ is very close to Is and close to p, but fer from Ip, whereas iP^ is 
plose to Ip but very far from Is and p. Out of the first 8 variables three are close to 
Ip alone, two are close to both Is and p, one is close to Is alone, and two are close to p 
alone. It thus appears that Ip and Is do not contain the same information and that part of 
S the information contained in Is and in p seems to be common. 

5. 2 - Filtering of impedances P, S and of the density 

The three elastic properties can be reconstituted with a limited number of new 
seismic attributes. The number of synthetic variables used for reconstitution can be 
different from one set to the next : it depends on the inertia percentage to be 
10 reconstituted/filtered. 

5.3 - Supervised and non-supervised pattern recognition of the new attributes 

It is possible to interpret geologically the new seismic attributes through a facies 
analysis with a supervised or non-supervised pattem recognition algorithm. 

In supervised pattem recognition, it is necessary to select learning classes grouping 
IS together seismic events located in the vicinity of wells whose geologic interpretation 
has been carried out This methodology allows to classify all of the seismic events 
according to the geologic charact^s observed in the wells. 

In the present case, the geologic knowledge of the three available wells Wl, W2 
and W3, and the general geologic knowledge of the zone studied lead us to consider 
20 three classes : a first class corresponding to wells Wl and W2 represents turbidite sands, 
a second class corresponding to well W3 represents turbidite muds and a third class 
represents hemipelagites. Wells Wl and W2 showing close geologic characteristics, 
only the seismic events located aroimd well W2 are used, well Wl serving as a blind 
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test. The location of well W3 being uncertain, the seismic events used to characterize 
class 2 are located around a virtual well defined in the proximity of W3. Finally, since 
we have no well for class 3, a virtual well is defined jfrom our general geologic 
knowledge. Each class consists of 121 seismic events (11 x ll). These events are used 
5 to caUbrate a classification function which is then appUed to the whole of the seismic 
events. 

Figures 8A, 8B respectively show the assignments obtained and the most probable 
assignments, i.e. having an assignment probability greater, than or equal to 80 %. Well 
Wl, used in blind test, is correctly assigned to the same class as well W2. Well W3 is 
10 located at the boundary of classes 1 and 2. Class 2 is distributed along a channeling 
form. The map of Figure 8B shows that the results are stable : few events exhibit a bad 
assignment probability. 

In non-supervised pattem recognition, the learning classes are formed using the 
modes of a multivariate probability density function calculated firom the whole of tiie 

15 seismic events, in a space generated by the first eight new attributes. The seismic events 
of higher probability density that constitute the peaks of this density function are 
selected to carry out learning of the classification function. This function is then applied 
to the whole of the seismic events considered to assign them to one of the classes : 
Figures 9A, 9B show the assignments obtained, as well as the most probable 

20 assignments, i.e. having an assignment probability greater than or equal to 80 %. The 
maps obtained are very close to the maps of Figures 8A, 8B. Class 3 of the supervised 
analysis corresponds to classes 3, 4 and 5 of the non-supervised analysis. Wells Wl and 
W2 are assigned to the same class and well W3 is still located at the boundary between 
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classes 3, and especially 1 and 2. This convergence of results validates the interpretation 
of the supervised approach. 



